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Abstract 

This work is devoted to the numerical simulation of a Vlasov-Poisson model 
describing a charged particle beam under the action of a rapidly oscillating external 
electric field. We construct an Asymptotic Preserving numerical scheme for this 
kinetic equation in the highly oscillatory limit. This scheme enables to simulate the 
problem without using any time step refinement technique. Moreover, since our 
numerical method is not based on the derivation of the simulation of asymptotic 
models, it works in the regime where the solution does not oscillate rapidly, and 
in the highly oscillatory regime as well. Our method is based on a "double-scale" 
reformulation of the initial equation, with the introduction of an additional periodic 
variable. 

1 Introduction 

In this article, we are interested in the construction of numerical schemes for collisionless 
kinetic equations which involve rapid oscillations in time. Our study is done in the 
framework of a specific physical application, the case of a charged particle beam in 
the paraxial approximation, but our strategy can be applied to other highly oscillatory 
kinetic models, for instance in the physics of magnetized plasmas [15, 13, 14, 4, 5] for 
the guiding-center limit or the finite Larmor radius limit. 

Let us first present our model. The paraxial approximation of the Vlasov-Maxwell 
equations concerns stationary, non collisional, charged particle beams which display a 
predominant length scale, called the longitudinal direction, such that the transverse 
width of the beam is very small compared to the typical longitudinal length. The 
paraxial model is obtained by expanding the Vlasov-Maxwell model with respect to the 
ratio e > between the characteristic lengths in the transverse and in the longitudinal 
directions, we refer to [9, 10] for a derivation of this model. Here, following [3, 12, 22], 
we consider the simpler case of an axisymmetric beam (with zero angular momentum). 
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The paraxial Vlasov-Poisson model takes then the following form, in dimensionless 
variables, 

dtf + ^drf + {Efe + E^pp)d,r = 0, (1.1) 

where f^{t,r,v) is the distribution function of the particles, t > corresponds to the 
longitudinal position coordinate (the direction of propagation of the beam, denoted as a 
time here), r G IR+ is the radial component of the position in the transverse plane, and 
u € M is the radial velocity in this plane. The total electric field has two contributions, 
the self-consistent electric field Eje = Ef£(t,r) satisfying the Poisson equation in the 
transverse plane, written in cylindrical symmetry as 

-drirEf.) = [ rdv (1.2) 
r Jr 

and an applied electric field i^app, chosen as in [12] under the following form 

^app(t,r) = -J+a(^^)r, (1.3) 

where a is a given 27r-periodic function (the so-called tension function). This system 
is initially defined for r > but can be extended to r S M by using the conventions 
f%t, -r, -v) = f%t, r, v) and E{t, -r) = -E{t, r). 

To summarize, in this paper we consider the following one-dimensional Vlasov- 
Poisson system satisfied by f'^{t,r,v), where r G M and u G M, 

dtf + ^drf+(^Ef.-^-+a(^^r^d,f = 0, f{t = 0,r,v)=Ur,v), (1.4) 

Efe{t,r) = -[ sp%t,s)ds with p''{t,r) = [ f{t,r,v)dv. (1.5) 
f Jo Jr 

The initial data /o is a given smooth function. When there is no confusion we shall 
omit the subscript e to ease notations. 

The main purpose of this work is the construction of robust numerical methods for 
stiff transport equations of type (1.4) in the limit e — )■ 0. We seek a method that is 
able to capture the properties of the various scales in the considered system, while the 
numerical parameters may be kept independent of the stiffness degree of these scales. 
Contrary to collisional kinetic equations in hydrodynamic or diffusion asymptotics, 
collisionless equations like (1.4) involve time oscillations. In this context, the notion 
of two-scale convergence [1, 12, 23, 8] is well-adapted in order to derive asymptotic 
models. However, these asymptotic models are valid only when e is small. In this 
paper, we develop numerical schemes that are able to deal with a wide range of values 
for £. We construct a numerical method in the so-called Asymptotic Preserving (AP) 
class [17]: such schemes are consistent with the kinetic model for all positive value of 
£, and degenerate into consistent schemes with the asymptotic model when e — )• 0. 
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To do this, let us first put the stiff equation (1.4) into a filtered form by rewriting it 
in the adapted rotating frame. The characteristic equations associated with (1.4) read 

Jt{ v) ^ e"^ { v) ^ { Ef{t, r)+a (t/e) r 



where the matrix J is defined by 

J 



1 
-1 



Hence, introducing the oscillatory variable ^ G defined by 

^1 ^ _ ^-Jt/e ( r \ _ ( cos{t/e) -sin(t/e) \ ( r \ 
i2 ) W y V sm(t/e) cos(t/e) )\v ^'"^^ 

the associated filtered distribution function /^(i, '^i, ^^2) = f^it-iT^v) satisfies 

dj^{t,o + (^p(t,t/e,e) +^app(t/e,o) • (t,o = 0, ~f{t = 0,-) = /o, (1.7) 

where the vector field is the sum of the applied field 

' — ' / — Sill T \ 

£^app(r,e)=o(r)(eicosr + 6sinr) (1.8) 



cos r 



and of the self-consistent field defined by 

/ -sinr \ 1 f+°° ~ 

Ef(t,T,S.)=[ ] — -7 / / s f (t, s cos r — f sin r, s sin r + -y cos r) dsfii) 

^ V COST J r{T,0 Jo J^oo 

with r(T, C) = Ci cos r + ^2 sin r. 

Let us briefly describe the strategy we propose to deal with equations like (1.7). As 
a matter of fact, we embed the function f^{t, ^) into the family of solutions F^{t, r, S^) of 
an "augmented" kinetic equation, where we separate the two scales t/e and t. Assume 
indeed that solves the equation 

dtF' + (^^.(t,T,0 +^app(T,e)) • V^F^ = -^drF', (1.9) 

and that, additionnally, we have 

veeK^, F^(o,o,o = /o(e), (1.10) 

then it is readily seen that F^{t, t/e, satisfies the initial-value problem (1.7), so we re- 
cover f^{t, ^) = F^(t, t/e, ^). The point is, in this double-scale formulation (1.9) of (1.7), 
the stiffness is confined in the sole term —^drF^ in the right-hand side. Reinterpreting 
this singularly perturbed term as a "collision" operator in this collisionless context, we 
can obtain the asymptotic behavior of F^ (then of /^) by a Chapman-Enskog expan- 
sion. In turn, this suggests a systematic method to construct Asymptotic Preserving 
numerical schemes, based on a micro-macro decomposition of F, see [20, 2, 7]. 
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This paper is organized as follows. In Section 2, we present the double-scale for- 
mulation in a general framework, and perform in subsection 2.1 the Chapman-Enskog 
expansion of F. We then discuss in subsection 2.2 the crucial question of the choice 
of the initial data i^^(0, r, ^) for this augmented kinetic equation (1.9). In this subsec- 
tion, we state the main (formal) theoretical result of this paper, in Proposition 2.1. In 
subsection 2.3, we compute explicitely the averaged equations for our problem in the 
linear setting (when self-consistent interactions are neglected). Then, in Section 3, we 
present our AP numerical scheme. In subsection 3.1, we introduce the scheme, which 
is a second order (in time and space) Eulerian numerical scheme. In subsection 3.2, we 
prove formally that this scheme is Asymptotic Preserving at the limit e ^ 0. In sub- 
section 3.3, we show how the micro-macro decomposition method enables to construct 
AP schemes in more complicated situations, such as the diffusion limit. Finally, the 
last Section 4 is devoted to a series of numerical tests which characterize the properties 
of our scheme. 

2 Double-scale formulation of the oscillatory equation 

In this section, we introduce a general strategy in order to deal with highly oscillatory 
problems under the form 

dtf' + A{t,t/e,th = 0, f(t = 0,-)=/o, (2.1) 

where the unknown is the distribution function (t,^) G x M'^ i— )• f'^{t,^) G M and 
the vector-field (t, t, ^,/) i— )• A(t,r, ^, /) G M is a functional which is P-periodic with 
respect to the variable r G T (T denotes the torus W/PZ). Our target equation (1.7) 
is under the form (2.1), with d = 2 and 

A{t, r, f) = {Ef{t, T, + ^app(r, O) ■ V^/. 

We now introduce the following "double-scale formulation" 

dtF' + A{t,T,C,F') = -^drF', (2.2) 

where the unknown is the function [t, r, ^) G M+ x T x R"' i— F^{t, r, ^). This problem is 
an augmented version of (2.1). Indeed if a function F^{t,T,S^) solves (2.2) and satisfies 
additionally 

VCGM^ F^(0,0,0 = /o(0, (2.3) 

then by differentiating F^{t,t/e,(,) we obtain that f^{t,(,) := -F^(t, i/e,^) satisfies the 
initial-value problem (2.1). 

It is important to note that (2.2), (2.3) is not sufficient to uniquely determine the 
function F^. Indeed, (2.3) is not a Cauchy condition for (2.2). The question of choosing 
a "good" initial condition F(0,r, ^) = Fq{t,S^) for all (r, ^) G T x M"^ is a delicate issue 
and is discussed in subsection 2.2. In fact, we will see - in a formal setting - that there 
is a unique way (up to order 0{e'^) terms) to define Fq in order to get a smooth function 

(t,r,e,e) G [0,t/w] X T X M-^ X [0,eo[^ F^t,T,0 
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that satisfies (2.2) and (2.3). Here tfinai > is a fixed final observation time and eo > 
is arbitrary. The important point here is the assumed regularity with respect to e when 
this parameter goes to zero. 

More precisely, our aim is to ensure that the function and its two first derivatives 
dtF^ and d'^F'^ are bounded. Roughly speaking, this regularity is a constraint that 
prevents a dependency of F'^ in the fast variable t/e (up to order O(e^) terms), and F^ 
will "only depend" on t, r and Under this condition, one can pretend that we have 
succeeded in separating (up to order O(e^) terms) the two scales t and t = t/e that 
were initially in (2.1). The main result of this section is Proposition 2.1. 

2.1 Chapman-Enskog expansion 

In this subsection, we analyze formally the behavior of (2.2) when e — t- 0, assuming 
that its solution F^ is smooth enough. To this aim, we carry out the Chapman- 
Enskog expansion of this function. Consider the following linear operator, defined for 
all periodic (regular) function r E T i— )• /i(r) by 

Lh = drh. 

This operator is skew-adjoint with respect to the L^(T) scalar product and (2.2) can 
be rewritten 

dtF' + A{t,T,i,F') = --^LF'. (2.4) 

The kernel of L is the set of constant functions and the L? projector on this kernel is 
the average 




where |T| = P is the measure of T. 

Moreover, L is invertible in the set of functions with zero average and, if Jj, h{T)dT = 
0, we have 

(L~^/i)(t) = (I-n) f h{a)da= [ /i(a)dcj + ^ [ ah{a)da. 

Jo Jo FI Jt 

Performing the Chapman-Enskog expansion of F^{t,T,^) consists in writing 

F%t,T,C) = G%t,C) + h%t,r,0 with G'{t,0=Il{F%t,T,0) (2.5) 

and deriving asymptotic equations for and /i^ when e — t- 0. As we said, we proceed 
at a formal level, and the rule that we follow in this analysis is that F^ is assumed to 
be smooth with respect to all its variables (in particular with respect to the parameter 
e which can be very small). 

Inserting the decomposition (2.5) into (2.4) leads to 

dtG" + dth^ + A{t,T,tG' + h") = --Lh'. (2.6) 
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Averaging this last equation with respect to r (i.e. applying 11) yields, since Hh^ = 0, 



dtG' +n{A{t,T,^,G' + h')) = 0. (2.7) 
Then, from (2.6) and (2.7) we deduce that satisfies 

dth' + (I - n) {A{t, r, ^, + h')) = -\Lh'- (2.8) 
Now, from (2.8) and the fact that belongs to the range of L, we deduce that 

= -eL-i {dth' + (I - n) {A{t, r, e, G' + h'))) . (2.9) 
Hence, using our smoothness assumption and in particular that we have dtF^ = 0(1), 



d^F'^ = 0(1) (hence and /i^ have also bounded derivatives), we deduce from (2.9) 
that 

= 0{e) and dth = 0(e). 
Prom these estimates and (2.7), we deduce a first approximate equation satisfied by 

dtG' + I[A{t,T,i,G') = 0{e). (2.10) 

Next, using again (2.9), we obtain an expression of in terms of G^ , up to a small 
remainder: 

= -eL-i(I - Tl)A{t, T, i, G") + 0{e^) (2.11) 

and this expression, together with (2.7), enables to derive the following equation satis- 
fied by G^ up to second order terms: 

dtG' + TlA{t, r, G') - eli (dfA{t, r, G') (I - Ii)A{t, r, G'))) = 0{e^). 

(2.12) 

Finally, the function can be deduced from G^, up to second order terms, by using 
(2.5) and (2.11): 

= G^ -eL-i(I-n)^(t,r,^,G^) + 0(e2). (2.13) 
2.2 Discussion on the initial data and main result 

In the previous subsection, the Chapman-Enskog expansion was performed formally 
under a regularity assumption on F'^ . In this subsection, we reverse the argument and 
deduce from these expansions a Cauchy data for (2.2) that ensures that F^ is regular 
enough (up to order 0{e^) terms). 

A natural initial condition for (2.2) can be deduced from (2.13). Indeed, by evalu- 
ating (2.13) at t = 0, one gets 

F^(0, r, = G^(0, - e{I - H) [{I - n)^(0, s, G^(0, i))ds + 0{e^) (2.14) 

JO 

and then, by taking this equation at r = and by using (2.3), 

/o(e) = G^(0, i) + en [{I - n)^(0, s, e, G^(0, i))ds + 0{e^). (2.15) 



6 



By substracting these two identities (2.14) and (2.15), one gets 

F'{0, r, = /o(0 n)A(0, s, G'{0, 0)ds + 0{e^). (2.16) 

JO 

Moreover, from (2.15), one deduces G^(0, ^) = /o(C)+C'(e), which can finally be inserted 
into (2.16) and yields 

i^^(o, T, e) = /o(e) - e f\i - n)A(o, s, e, um^ + (2.17) 

The correction term in e is important here and, as we show further, will guarantee 
that F^{t,T,^) does not oscillate in time. By analogy with boundary value problems in 
collisional kinetic theory (see [27]), one can interpret this term as "boundary corrector" 
(where the boundary is the initial time t = 0). The interesting point in our case is that 
we do not have to assume that the initial data is well-prepared since, as we said in the 
introduction of this section, we have a degree of freedom on Fq which is not totally 
prescribed. We have then the possibility to enforce that (2.17) is satisfied (see (2.19)). 

Let us formulate in the following proposition the main result of this section. 

Proposition 2.1 (formal) Let F^{t,T,S^) be the unique solution of (2.2) subject to the 
initial condition 

V(r, e) G T X F^(0, r, = Foi^, (2-18) 

with Fq defined for all t ^ T and ^ £M? by 

FSir,0 = fo{0-e [\l-n)AiO,s,tfo{0)ds. (2.19) 
Jo 

Then we have 

F%t, T, = G%t, - eL-\l - U)A{t, r, t &{t, 0) + 0{e^), (2.20) 
where G^(t,^) is the solution of the initial-value problem 

dt& + IiA{t, r, &) - eli {dfA{t, r, &) - I[)A{t, r, i, = 0, (2.21) 

G^'(0,0=nFo^(r,0 = /o(0-^n r{I-U)A{0,s,^,foims. (2.22) 

Jo 

Remark 2.2 Since (2.18) and (2.19) imply (2.3), one can recover the solution to 
the oscillatory equation (2.1) by setting 

f{t,0 = F'{t,t/e,0. 

Moreover, whereas the term of order e in f varies rapidly in time (it depends on t and 
t/e), the corresponding term in the function F^ is smooth - since t replaces the variable 
t/e- and then easier to compute numerically. Our Asymptotic Preserving numerical 
method is constructed on the double- scale formulation (2.2) instead of {2.1). 
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Proof of Proposition 2.1. Let F" be the solution of (2.2), (2.18) and let be the 
solution of (2.21), (2.22). Denote 

with 

If = -eL~^{l - Ti)A{t, r, i, G^) 

and where is a bounded corrector that is defined below (see (2.24)). Proving the 
Proposition amounts to proving that 

F'{t,T,i)-F'{t,T,i) = 0{e''). 

By substracting (2.18) and (2.22), one gets 

F^(0,r,0-G^(0,0 = -e(I-n) r(/-n)yl(0,s,e,/o(e))rf5 

JO 

= -e(I-n) r {I - Ii)A{^,s,i,G'{^,i))ds + 0{e'^) 
Jo 

= -eL~\l - U)A{0, s, ^, G^(0, 0)ds + 0{e^) 

This gives 

F^(0, r, - F%0, r, = -e\'{0, r, + 0{e^) = 0{e^). (2.23) 

Let us now derive an approximate equation satisfied by F^{t,T,£,). By inserting F^ in 
the equation (2.4), one gets 

dtF' + ^LF' + A{t,T,tF') 

= dt& - eL-\l - U)dfA{t, r, ^, G^)dt& - eL-^{l - Il)dtA{t, r, &) 
-{\-Yi)A{t,T,i,&)+eLx' 

+A{t, T, e, G^) - edfA{t, T, e, G^) - n)^(t, r, &)) + 0{e^) 

= -IlA{t, r, e, G^) + en {dfA{t, r, G^) (l-1(I - li) A{t, r, G^)) ) 

-eL-i(I - n)5/^(t, r, G^)nA(t, r, ^, G^) - eL-^{l - Ii)dtA{t, r, G^) 
-(I-n)A(t,r,e,G^)+eLx^ 

+^(t, T, e, G^) - eayA(t, r, ^, G^) (l-1(I - n)^(t, r, G^) + ©(e^) 
= eLx' - e(I - n) {d}A{t, r, e, G^) - n)A(t, r, G^) ) 

-eL-i(I-n)5M(t,r,e,G^) 

-eL-i(I - n)5/^(t, r, G")nyl(t, r, + O(e^) 
where we used (2.21) in the second equality. Hence, by defining the corrector x^ 

X' = L'^[{\-Ti){dfA{t,T,i,&){L'\l-Yi)A{t,T,i,& 

+L-i(I - Ii)dtA{t, r, G') + L~\l - U)dfA{t, r, G')UA{t, r, (2.24) 
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one gets finally 

dtF' + ^drF' + A{t,T,tF') = 0{e^). (2.25) 

We shall now conclude by integrating the characteristics associated to this equation. 
Let 

w{t,r,0 = {F'-F'){t,r + t/e,0. 

Prom (2.2) and (2.25), one deduces 

dtw = -A{t, T + t/e, F%t, T + t/e, C)) + A{t, t + t/e, F^{t, t + t/e, 0) + 0{£^). 

Hence, using the estimate (2.23) at the initial time, a Gronwall lemma yields formally 
(recall that A is periodic with respect to r) 

w{t,r, = 0(8^) 

for t e [0,r], r G T, ^ G and, finally, one has proved that 

(f^ - F^) {t, T, = w{t, r - t/e, = ©(e^). 

The proof of Proposition 2.1 is complete. ■ 

Remark 2.3 In fact, this averaging procedure can be pushed forward to higher orders 
in £ by iterating further the Chapman- Enskog procedure. Other approaches may he 
used to obtain formally higher order averaged models for F^ , under a higher order 
initial condition, see for instance the approach developed in [24, 26] which is widely 
used in the context of ODEs. However, the purpose of this paper being to build an AP 
numerical method for our problem, we stop this construction at order O(e^). We also 
refer to [6] for a presentation of the so-called stroboscopic averaging in a way which 
is very close to the method introduced here. Indeed, in [6], a systematic construction 
of high order averaged models for oscillatory equations such as (1.7) is based on the 
transport equation (2.2). It is proved in this paper that, for any fixed integer N > 0, 
the solution of (1.7) can be written under the form (omitting the dependencies in ^ for 
simplicity and assuming that A{t,T,f) does not depend on t) 

f{t) = ^.^'^ {t/e, G^'^(t)) + 0(e^+^), (2.26) 

where G^''^ (t) satisfies an autonomous averaged equation of the form df G = j4|v^ (G) 
with G^'^ i/d) = /o and where {t, f) i— )• $''''^(r, /) is a close-to-identity mapping which 
is 2TT-periodic with respect to r and satisfies $^'^(0,/) = /. The link with our con- 
struction is the following. If we choose Fq{t) = ^"^'^ (r, /o) as initial data for (2.2), 
then the stroboscopic averaging result says that F{T,t) = (r, G'^'^(i)) + 0{£^~^^), 
i.e. F{T,t) is smooth, up to 0{£^~^^) terms. This gives the natural generalization of 
our initial data (2.19) in order to get higher order estimates. Of course, one can check 
that, for N = 1, 

f (r, /o) = fo-e Hi - U)A{s, /o(0)ds. 

^0 
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2.3 The case of a linear transport equation 

In this subsection, we compute explicitely the initial condition Fq and the averaged 
system in the special situation of the following linear transport equation in dimension 
d = 2: 

dtf' + E{t, ■ Vgf = 0, fit = 0, •) = /o, (2.27) 

where the field E{t,S^) = ^ ^ ^® given and divergence-free. This equation is 

under the form (2.1) with 

A{T,^,f) = E{T,0-yd- (2-28) 
In particular, when the self-consistent Poisson field Ef is neglected, the filtered equation 
(1.7) associated to the paraxial beam model (1.4) is under this form, with E[t,^) = 
Ea.ppiT,C) defined by (1.8) (it is a divergence-free vector field). 

In this linear case, the following proposition is a variant of Proposition 2.1. 

Proposition 2.4 (formal) Assume that A takes the form (2.28). Let F^(t, r, ^) be 
the unique solution of (2.2) subject to the initial condition F^(0, r, ^) = Fo(t, ^) with 
Fq defined for all t £ T and ^ G 5y 

Fo(r, = /o - e j\l - n)i^(s, Ods^ ■ (2-29) 

Then we have 

F'it, r, = & {t, C - eL-' (I - n)i?(r, 0) + O(e'), (2.30) 
where G^(t,^) is the solution of the averaged transport equation 

dt& + + eS(^)) • = 0, (2.31) 

= fQ(i - eU j\l -Ii)E{s,i)ds^ , (2.32) 

and where E^^^ = HE and E^^^ = J~^V^P is the vector-field associated with the 
Hamiltonian 

-Dii) = ^ / [(^ - n)^2] (r, j\l - ^)Ei {s, OdsdT. 

Remark 2.5 This result is the Eulerian version of an averaging theorem formulated 
directly in terms of the characteristics equations associated to the vector field E{t,^). 
Indeed, consider the flow H associated to the averaged vector field: H(i,to)Co) solves 

^ = (H) + (H), H(to, to, Co) = Co. 

Then we have F^(t,r, C) = /o(H(t,r, C)) + 0{e^), where H is defined by 

E{t,T,0 = (^l-enJ\l-U)E{s,-)ds^ (e (^0,t,( - e j\l - U)E{s,Ods^'j . 
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Remark 2.6 The averaged equation (up to the the order 0{e'^)) shares the geometric 
structure of the initial equation (2.27). Indeed, since E is divergence-free, so is 
and if E is Hamiltonian, with Hamiltonian H{t,^), then E^^^ is Hamiltonian, with 
Hamiltonian given by H^^^ = IIH . Moreover, the correction eE^^^ is always divergence- 
free and Hamiltonian. 

Proof of Proposition 2.4. The initial data (2.29) and (2.32) can be deduced from 
(2.19) and (2.22) by a Taylor expansion, up to order O(e^) terms: one has indeed 

/o (e + eB{T, 0) = /o(0 + eB{T, ■ V^foiO + 0{e^). 

Similarly, the change of variable (2.30) can be deduced from (2.20) by a Taylor expan- 
sion. Moreover, we have clearly IlA{t,T,^,G^) = E^^^ ■ V^G^. Hence, to end the proof 
of the proposition, we simply have to compute the first order correction in the equation 
of given by Proposition 2.1, i.e. the operator 

G ^ -U{dfA{t,T,C,G){L-\l-U)A{t,r,tG))) 



~ [ E-V^ {L-^il - U)E ■ V^G) dr = • (BV^G) 



where we used that E is divergence- free and where B is the 2x2 "diffusion" matrix of 
components 

' E,L-^[iI -U)E,]dT, i,j = 1,2. 



1 

W\ 



In fact, this matrix 
we have 



inherits the skew-symmetry property of L. Indeed, for all i,j, 

= -^J'^rnil -U)L-'[{I -U)Ej]dT 
= -^J'^[{I -U)E,]L-'[iI -U)E,]dT 
= [ L-\I - Il)Ei\ {I - U)EjdT = -B,-,. 

I ^1 JT 

Hence, setting D = Bi^2 = —^2,1, the "diffusion" term • (BV^G) can be simplified 
as 

• (BV^G) = %(P%G) - %(P%G) = (%P)%G - (%P)a5,G, 

which is the desired result. Note that the first order model is a pure transport equation 
and does not include second order derivative. ■ 



Explicit calculations in our example. 

Let us compute explicitely the approximate model in terms of the Fourier coefficients 
of E. From now, the period is taken as |T| = P = 27r. Introduce the decomposition of 
the two (real-valued) components of the vector field E on the Fourier basis: 

E,{T,C) = Y,^kAOe"'^ for J = 1,2, 
fcez 
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with A_kj = Afcj for all € Z and j = 1,2. Then direct calculations yield 

Ef=Aoj for j = 1,2 and V = 2Im l^k^iA^ . (2.33) 

We now calculate the quantities defined in Proposition 2.4 in a specific example that 
we use later for numerical experiments. In the beam model (1.7), if we neglect the 
Poisson field, then we have E = -Eapp defined by (1.8). Choosing a(r) = cos^(2r), one 
computes from (1.8) the Fourier coefficients of Ei and E2: 

El = ^(-46 + (36 + ^6)e''"-26e''" + (6 + ^6)e''" + c.c.) 
= ^ (46 + (36 - i6)e'*" + 26 e^^^ + (6 - i6)e'^*" + c.c.) . 
Hence, we obtain by simple integrations 

u[ (l-U)E = Do^ and L~^{1 - U)E = Di{t)C, (2.34) 
Jo 

with 

1 / -1 



12 V 1 

_ J_ / 3cos(2t) + cos(6r) 9sin(2r) - 3sin(4r) + sin(6T) 

^ ~ 48 V 9sin(2r) + 3sin(4r) + sin(6T) -3cos(2r) - cos(6t) 

and also, from (2.33), we obtain that the averaged vector field (up to order O(e^) terms) 
is the following Hamiltonian vector field: 

with 

The averaged equation (2.31), (2.32) for is thus the equation of a rotation in the 
phase space and has an explicit solution: 

G%t, = G%0, e'^'O = fo ((I - eDo)e'^-'^) . (2.36) 

We have thus analytic expressions for the solution of the limit model as e — t- and 
also for the solution of a next order approximation, which are then easy to implement 
numerically. The solution of the limit model reads (see [12]) 

i^iimit(t,r,0 = /o(e*""-^0 (2.37) 
and the solution of the second order model will be 

i^sccondordcr(t, T, = fo ((I " e^o)e*"-^(I " eDi{T))C) . (2.38) 
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This last relation is obtained by using successively (2.30), (2.36) and recalling that Dq 
and Di are given by (2.34). Indeed, one has 



= G%0,e^'^-^{I -eDi)C) + 0{e^) 

= fo ((I - eDo)e''''{l - eD,{T))C) + 0{e^). 

3 Asymptotic Preserving numerical schemes 

In this section, we construct some Asymptotic Preserving numerical schemes for (2.2), 
hence for the original problem (1.4). Let us insist on the fact that we do not base 
the construction of our numerical method on the approximate models derived in the 
previous section, since we want a method which is efficient for the regimes where e 
small and where e = 0(1). 

Recall that, in order to solve the filtered equation (1.7), we have introduced the 
augmented equation 

dtF' + E{t,T,0-y!;F' = -\drF', (3.1) 
where we denote for simplicity the field (which depends on the unknown F'^) by 

s(t,T,e) = ^F^^(t,r,e) + ^app(r,e). 

After the asymptotic analysis in the previous section, and according to Proposition 2.1 
(see also Proposition 2.4), we know (see (2.29)) that a suitable initial condition for this 
problem is F(0,r, ^) = Fo(r, ^) with 

i^o(r, i) = h{i-e j\l - H)E{Q, s, C)ds^ . (3.2) 

Note that this choice is asymptotically close to (2.19), up to order 0{e'^) terms, but is 
preferable since it garantees the positivity of the initial distribution function. Under 
this choice, we know two important facts: 

- one recovers the solution of (1.7) by f{t,(,) = F^{t,t/£,^), 

- the function is smooth and, up to terms of order O(e^), does not oscillate. In 
particular, its derivatives dtF^ and d'^F^ are bounded when e — t- 0. 

In order to emphasize the role of the choice of the initial condition Fq , in our numerical 
experiments we will also test the most simple choice: 

Fo{r,0 = fo{0- (3.3) 

This choice only garantees that F'^ = + 0{e): we show below that, with this initial 
data, the numerical method will capture the right limit, but not the details of order 
0{e). In the sequel, the initial condition (3.2) is referred to as "with correction", and 
the initial condition (3.3) is referred to as "without correction". 
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3.1 The numerical scheme 



In this subsection, we present our AP numerical scheme. Due to the lack of relaxation 
or diffusion operator, very fine structures and filamentations can be observed which 
forbid the use of low order numerical methods. Hence a second order (in time t and 
phase space ^) finite difference discretization is applied to (3.1), which is based on a 
Lax-Wendroff-Richtmyer numerical scheme (see [25, 16]). 

First, we introduce the time discretization = nAt with n G N and the time step 
At. The phase space discretization is uniform so that the domain [— ^max> ^max]^ is 
meshed by = -^max + and ^2,i = -'^max + jA^ for f, j = 0, . . . , iV - 1 and 

= 2^max/-^) ^ being the number of points per direction. For the direction r, we 
also use a uniform mesh of size At, so that T£ = iAr, for i = 0, . . . , A^r — 1? Ar = 
211 /Nt-. Denoting = (^i,j,^2j)i the discrete unknown is then F^j ^ ~ F^{t"',Ti,^ij). 
In the following description, we keep the r variable continuous in order to focus on 
the discretization in the and ^2 directions. In practice, since periodic boundary 
conditions are considered in this direction r, the fast Fourier transform is very efficient 
for this variable. At the boundary of the phase space domain in ^, zero inflow boundary 
conditions are prescribed. 

We then introduce the flux in ^ which approximates (E'" • 'V^)F^j by centered finite 
differences: 

rpn Tpn rpn rpn rpn rpn rpn jpn 

, ^2,i,j+l^i,j+l ^2,i,j-l^i,j-l 



2A^ 2A^ 
and we also consider the following four-points average 

F^^^ = [Fix^, + + + /4. 

A first step on At/2 is performed to get intermediate unknowns F^^^^'^ 



pn+1/2 pn _At At n+1/2 



The second step reads 

pn+l ^ pn _ ^^^n+l/2^^„+i/2) _ ^Q^^F,^^ + ^.-f 1). (3.5) 



Standard results (see [16, 25]) say that this numerical scheme is second order in time 
and phase space for all fixed e > 0. 

Recall now that the model is nonlinear due to the presence of the self-consistent 
electric field Ep. Let us explain how we update the field Ep'^^, once is known. 

The inversion of the Poisson equation is easier in the original variables (r, v) than in 
the variables since it takes the simple form (1.5) of an ODE in the r variable. At the 
continuous level, coming back to (r, v) can be done easily by introducing the function 

f{t,T,r,v) = F(t,r,6,6), with ( f]= e""-^ ( . (3.6) 



6 7 V ^ 
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It is not that simple at the discrete level. Indeed, if {ri,Vj) is the mesh in {r,v), then, 

for all given r^, the points e~'^^'^ ^ ^ ^ do not necessarily coincide with mesh points 

^ij. To evaluate f"'~^^{Ti,ri,Vj), we thus need an interpolation algorithm in dimension 
2. Since this interpolation is done at each step, we choose a simple linear interpolation 
algorithm. Then, once we known the values fi'j'^ for each Te, it is easy to compute 
the Poisson field EJ~^^ by integrating (1.5). To deduce Ep'^^ on the ^ mesh, another 
interpolation step is required. Finally, we also remark that our algorithm in two steps 
imposes to predict the advection field E at time t*^"*"^/^, so a Poisson field evaluation is 
needed also before computing the flux 

At the final time t/maZ of the simulation, we come back to the solution of our initial 
problem (1.4) by setting f {t final, r,v) = F{tfinaht final 1^,0, so a last interpolation 
algorithm in the two-dimensional (r, v) variable is needed, as well as in the r variable 
(since t final does not necessarily coincide with a discrete Ti). 

3.2 Asymptotic Preserving property 

In this subsection, we check formally that the numerical scheme presented above is 
Asymptotic Preserving, as announced. Thanks to the implicitation of the stiff term jdr, 
the only stability condition will be a standard CFL condition of the form At < CA^. 
In the sequel, we consider for simplicity that At ~ A^. We have already seen that, for 
fixed e > 0, this scheme is consistent (and of order 2) with the equation (3.1). We now 
have to examinate its behavior when e — )■ 0. 

It is convenient to analyse the asymptotics of numerical schemes written with the 
micro-macro decomposition technique, which was developed in [20, 2] as a flexible 
method in order to construct Asymptotic Preserving numerical schemes for collisional 
kinetic equations. Remark that, here, we have rewritten (1.7) under the "collisional 
form" (1.9) (the operator dr plays the role of the collision operator). The micro-macro 
method consists in mimicking the Chapman-Enskog expansion and decomposing the 
unknown F^ into a macro part = YIF^ and the remainding micro part = (I— 11)^^. 
This micro part is small when e is small (but plays an important role when e is not 
small, ensuring the AP property). In fact, our scheme (3.4), (3.5) is already under a 
"micro-macro" form, thanks to the simple form of the operator L = dr- Indeed, it 
suffices to set 

= ^E^j, hi j = (I - n)Fjj 
to realize that our scheme is reformulated as follows: 

(3.7) 

hil'/' = - f (/ - u)<^i^{G- + hn - ^drK^"' , 

'G^f = G^, - Atn$^+'/'(G"+i/2 + /i"+i/2). 

(3.8) 
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We will proceed by an induction argument. From our choice (3.2) of initial data, 
we deduce that 

G° = /i° = -eL~\l - n)$°j.(G°) + 0{e^). 

Now, assume that we have proved that 

G" = h"" = -eL-\l - U)^lj{G'') + 0(e2 + eAt). 

On the one side, the micro part of the first step (3.7) gives 

from which we deduce that 

h^;^'/' = -eL-\l-U)^l^{Gn + 0{e') 

= -eL"i(/-n)$"+^/^(G"+i/2) + 0(e2 + eAi), (3.9) 

since = + 0{At) and E"+i/2 = + ©(At). On the other side, the micro 

part of (3.8) leads to 



/i^+i = -2eL-i(/-n)$"+^/^(G"+i/2) _;jn^0^^2) 

= -2eL~\l - n)^'"+^/^(G"+i/2) + _ n)$^^.(G") + 0(e2) 

= -eL-\l -U)<^l+\G''+^) + 0{£'^ + eAt), (3.10) 

which ends the induction proof. 

Let us now focus on the AP property. The macro part of (3.7) gives 

= G^i-^n$^^.(G"-eL-^(/-n)$^^.(G™)) + 0(e2At + eAt2). (3.11) 
If we now insert (3.9) into the second equation of (3.8), we then obtain 
G^+i = Gl^ - Atn$^+'/' + /i"+V2^ 

= Gl^ - Atn$:;+^/' _ _ n)$^+^/2(G"+i/2)) + 0{e^At + sAt^). 

(3.12) 

Passing to the limit as e — )■ (for fixed A^, At) in (3.11), (3.12) yields 

'g::^/^ =G^-^^$^^.(G"), 
G7+1 = G", - Atn$rf f G"+i/2 

*ij '^iJ 
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which is a Lax-Wendroff-Richtmyer numerical discretization of the hmit equation 



dtG + UE ■ V^G = 0. 

This proves that our scheme is Asymptotic Preserving. Furthermore, we observe that 
when e is small but not zero, up to a 0(e^At + eAt^) remainder, the numerical scheme 
(3.11), (3.12) is nothing but a second order Lax-Wendroff-Richtmyer numerical dis- 
cretization for the approximate asymptotic equation (2.21) of . Hence, accumulating 
the errors will yield HG" — G^||oo < Ce^ + CeAt < Ge'^ + GAt'^ (here G denotes a generic 
constant independent of e. At and A^). 
We have then, for all n, 

= -eL-i(/-n)$^j.(G") + 0(e2 + At2) 
= -eL^\l -U)E -V^G' + 0{e'^ +At^), 

where we recall that we have assumed At ~ A^. Finally, in view of (2.20), we have 
jpn _ Qn _^ j^n _ jpe ^ 0(^2 _|_ /^f'iy £g^j.^ ^j^jg analysis concerns the asymptotics 

e — )■ 0. For a fixed e > 0, we already know that our scheme is of order two in time and 
space, which means that there exists a constant K^e) > only depending on e and not 
on At such that \\F^ — -^"^11 < i^(e)At^. These two behaviors can be summarized in the 
following estimate 

-i^^lloo < Cmin {K{£)At^,e^ + At^) . 

This means that our scheme is in fact a second order Asymptotic Preserving in the 
following sense: 

- for all fixed e, this scheme provides a second order approximation of the original 
equation (2.2); 

- when e — 7- 0, this scheme degenerates into a second order approximation of the 
system (2.20), (2.21), which itself approximates the original equation (2.2) up to 
0{e'^) terms. 

3.3 Extension to the diffusion limit 

The micro-macro decomposition is not only a tool to analyze the limit e — t- (as in [20, 
21]), it is also a practical method that allows to extend the construction of AP schemes 
to more complicated situations, see [2, 18, 7, 19] for instance for collisional kinetic 
problems. Let us briefly present another oscillatory example that will be developed in 
a future work. For simplicity, we present this example in the linear setting of subsection 
2.3. We still consider (2.27), (2.28) but assume now that the average of E[t,^) in r 
vanishes: HE = (this is the case for the paraxial beam model if the forcing term 
a(r) has no Fourier component in the frequencies 0, 2 or —2). Then, the limit field 
E^^"^ in (2.31) vanishes and it is convenient to rescale the time variable in order to get 
a non trivial model at the limit. This amounts to considering, from the beginning, the 
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so-called "diffusion scaling" of (1.4) (even if the final model here will not contain any 
second order derivative): 

dtf + ^d,r + (a (-] --^] d,r = 0. 



In this case, the associated equation in F takes the following form, where the variable 
T stands for t/e^: 

dtF' + -E{t,0 ■ V^F' = -^drF'. (3.13) 

Our micro-macro scheme for (3.13) will consist in decomposing the discrete unknown as 
F^j = Gfj+hfj, where the macro part Gfj = HF^" and the micro part hfj = (I — n)Fj" 
are calculated by 



(3.14) 



K:^' = Kj n)$,,(i(G"+i + G") + _ ^driKj + K:^'). 

(3.15) 

Let us briefly discuss the limit of this scheme as e — ?■ 0. Since, initially, one has 
= 0{e) (see the discussion in section 2.2), it is readily seen that our semi- implicit 
scheme will propagate this property. For all n, one has h"' = 0{e), so the flux terms in 
the equations for G in (3.14) and in (3.15) are not singular. The last equation implies 
that 

hl+^ = -eL~\l - n)$,j(G"+^) + Oie") 

if this property holds true at step n. Hence, since it is true at step n = 0, it holds true 
for all n. Consequently, one deduces successively from the three first equations of our 
scheme (3.14), (3.15) that 

G"f = G^ + ^nc^,,(L"i(/ - n)c|.(G")) + 0{e), (3.16) 

h-+^'^ = -eL-\l - n)$,,(G"+V2) + o{e^) 

and 

G^+i = Gl- + MYi^i^j{L~\l - I[)<^{G''+^'^)) + 0{e). (3.17) 

Finally, if we disgard the remainders 0{e), the limit scheme (3.16), (3.17) is a Lax- 
Wendroff-Richtmyer scheme for the limit equation for G: 

dtG -U{E- V^{L~^{I - U)E • V^G)) = 0. 

The scheme (3.14), (3.15) is thus Asymptotic Preserving in the diffusion limit. 
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4 Numerical results 



In this section, we present some numerical results for the paraxial beam model (1.4), 
(1.5) described in the introduction. In particular, our aim is to validate the Asymptotic 
Preserving property of our scheme. For all the simulations, the function a in the 
applied electric field -Eapp defined by (1.3) is chosen as a(T) = cos^(2r). In the first 
series of tests, in subsection 4.1, we solve the complete Vlasov-Poisson model. Then, 
in subsection 4.2, we restrict our study to the linear case when the Poisson field is 
set to zero, and where the asymptotic models (the limit model and its e-correction) 
are explicit and can be solved analytically, which provide some additional reference 
solutions for small e's. 

The initial condition for (1.4) is the same for all the simulations. It is taken as a 
Gaussian in velocity multiplied by a regularized step function in r: 



fo{r,v) 



V2 



--X[r) exp 



ira 



-^), X(r) = ^erf 



r + 1.2 
0.3 



—erf 
2 



1.2 



0.3 



(4.1) 



with a = 0.2. For all the simulations, the space- velocity domain is {r,v) E [—4,4]^. We 
represent on Figure 1 this initial data for {r,v) £ [—2,2]^. 





-fO(r,0) 



Figure 1: Plot of the initial data /q. Left: 2D plot of the function in the (r,v) space 
(zoomed for {r,v) G [—2,2]^). Right: the two curves r /o('^, 0) and v i— )■ fo(0,v). 



Let us list the numerical methods which are tested below: 

— our numerical scheme (3.4), (3.5) with the initial data Fq given by (3.2), contain- 
ing the 0{e) correction term, will be referred to as AP with correction; 

— the same numerical scheme (3.4), (3.5), but with the initial data (3.3), without 
the correction term, will be referred to as AP without correction; 

— a splitting method for the initial, non filtered equation, (1.4), (1.5): we apply a 
second order time-splitting method (Strang splitting) for (1.4), that we split into 

dtf + ^drf = and dtf + Ufe -!-+a(^\r) ^-f" = 0' 
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each split equation being solved by a spectral method based on fast Fourier trans- 
form; this method will be referred to as the splitting scheme; 

- in the linear case (see subsection 4.2), we have the analytic expression (2.37) for 
the exact solution -Fumit of the limit model as e — t- - referred to as limit model - 
and we have (2.38) for the solution Kjccondordcr of the limit model with the first 
correction in e - referred to as second order model. 

For all the simulations, the number of discretization points in the r direction is = 64, 
hence the derivative dr and the integrals Jq are calculated with a spectral accuracy. 
The strategy for the choice of the time step is the following. For the two AP schemes, 
the time step is taken independently of e, it only has to satisfy the stability CFL 
condition related to our Lax-Wendroff-Richtmyer scheme, i.e. we always choose At = 

max ni'iX I ii/ 1 , with ^max — 

4 and A^ = 2^ 

max/-^5 ^ being the number of points 
in the (or in the ^2) direction. For the splitting scheme, we have to adapt At 
proportionally to e. The limit model and second order model are analytic and do not 
require any time discretization. 

4.1 The Vlasov-Poisson model for the beam 

Our first series of simulations concern the full model (1.4), (1.5) or its filtered equivalent 
version (1.7). 

Qualitative results for different regimes in e: Figure 2 

Let us start with a few qualitative results. We first show some 2D plots of the function 
at the same final time tfinai = '?r/4, for the three values e = 1, e = 0.25 and e = 0.01. 
We compare in Figure 2 the numerical solution obtained by AP with correction (here 
N = 128), to the reference solution computed with the splitting scheme with an adapted 
small time step. The time step for our AP scheme is At = 0.02 for the three values of 
£. These plots show a good agreement between our solution and the reference solution: 
the scheme AP with correction is able to capture all the regimes in e. 

Long time behavior and filamentation: Figure 3 

Now, we show that our AP scheme is able to capture very thin structures, with a 
numerical cost independent of e. On Figure 3 we plot the numerical solution obtained 
with the scheme AP with correction (with N = 512), for a very small £ = 0.001 and 
for different times t = tt, t = Att, t = 7tt and t = lOvr. We observe the filamentation 
due to the self-consistent Poisson field effect (compare to Figure 11 below, obtained at 
t = 2tt without the Poisson field). 

Numerical verification of the order 2 uniform accuracy with respect to e: 

Figures 4, 5 and 6 

Let us now proceed to more quantitative tests. We plot on the three next figures 
the relative L?' error between the numerical solutions computed with different schemes 
and a reference solution (computed with tiny time and space steps). The final time 
[t = 7r/16) is fixed. 
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e = 0.01, AP witi correction 
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e = 0.01, spUtting scheme 



Figure 2: 
t = 7r/4. 
with the 
£ = 0.01. 



2D plots for {r,v) G [—2,2]'^ of the numerical solutions f^{t,r,v) at time 
Left column: computed with AP with correction. Right column: computed 
spUtting scheme. Top line: e = 1. Middle line: e = 0.25. Bottom line: 
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For Figure 4, the solution is computed with the scheme AP with correction. On the 
left part, we represent (in logarithmic scales) the error as a function of the time step 
At, for different values of e (from e = 1 to e = 10~^): the slope is always close to 2 and 
the curves are very close together, indicating that the error is almost independent of £. 
This independence is confirmed on the right part of the figure, where we represent the 
error as a function of £, for different values of At: all the curves are nearly horizontal. 
These curves indicate that the error produced by the scheme AP with correction is of 
the form CAt^, with C independent of e. This proves experimentally the second order 
Asymptotic Preserving behavior of our scheme. 

For Figure 5, the same tests are done for the scheme AP without correction, i.e. 
for the scheme (3.4), (3.5) with the initial data Fq(t,(,) = /o(0- On the left part of 
the figure, we observe that the scheme behaves at an order 2 scheme for e = 0(1) 
(e = 1, 0.5 or 0.1) or for small values of e (less than 10"'^). But for intermediate 
regimes, the curves are more chaotic. On the right part of the figure, this feature is 
even more obvious: without the correction of the initial data, our scheme behaves well 
for £ = 0(1) and for e very small (in fact, when the observed error is greater than 
e), but not for intermediate regimes. This shows that this initial correction is really 
needed and this validates numerically the analysis done in Section 2. 

For Figure 6, the same tests are done for the splitting scheme (well resolved in 
space, we only observe the error in the time step). On the left part of the figure, we 
observe that, for all fixed e, the Strang splitting scheme is of order 2 but the important 
fact is that the error strongly depends on e: the smaller is e, the smaller must be the 
time step to maintain a constant error. We also observe this feature on the right part 
of the figure. Experimentally, one can estimate that the error for the sphtting scheme 
is of the form C{At/e)^. 

Evolution of an RMS quantity and observation of the oscillations in time: 

Figures 7, 8, 9 and 10 

Let us now observe the evolution in time of a Root Mean Square (RMS) quantity 
associated to the filtered distribution function f^{t,^i,^2)'- 



Note that, due to the filtering, this quantity does not oscillate at the limit e = 0, and 
only the corrective terms for e > are rapidly oscillating. On Figures 7, 8, 9 and 
10, we represent respectively, for e = 0.05, e = 0.025, e = 0.01 and e = 0.005, the 
time history of RMS{t) computed by the AP scheme with and without correction, and 
compare these numerical solutions to a reference solution. In all these simulations, 
we take = 128 and At = 0.02. In particular, for e = 0.01 and e = 0.005, the 
time oscillation is not resolved by this time step. However, in all the cases, one can 
observe that the solution obtained by the scheme AP with correction fits surprisingly 
well with the reference solution, see in particular the zooms on the right part of each 
figure: the red circles, which represent the only calculated points, are on the black 




(4.2) 
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Figure 4: Plot of the relative error for the scheme AP with correction. Left: error 
as a function of At for different e. Right: error as a function of e for different At. 
Conclusion: the scheme is of order 2 and the error (nearly) does not depend on e. 




Figure 5: Plot of the relative error for the scheme AP without correction. Left: 
error as a function of At for different e. Right: error as a function of e for different At. 
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Figure 6: Plot of the relative error for the splitting scheme. Left: error as a function 
of At for different e. Right: error as a function of e for different At. Conclusion: the 
error behaves like C(At/e)^. 



reference curves, even when the oscillation is not resolved. This is another proof of the 
Asymptotic Preserving property of our scheme. On the contrary, one observes that the 
solution obtained with the scheme AP without correction is less accurate: it converges 
to the right limit as e — ?• but it is not able to correctly give the details of order 0{e). 




Figure 7: Time history of RMS{t) for e = 0.05, computed by AP with correction and 
AP without correction. On the right: zoom of the left figure. 
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4.2 The linear case 

For all the numerical tests presented in this subsection, the self-consistent electric field 
is neglected in the model: we now set Ef = in (1.4). We are thus in the situation of 
the linear model of Section 2.3, for which we have analytic expressions for the solution 
-^limit of the limit model and for the solution Fgecondorder of the second order model, 
respectively given by (2.37) and (2.38). 

Qualitative results for two regimes of e: Figure 11 

As above for the nonlinear model, let us start with a few qualitative results. We first 
represent the 2D plot of the solution of the linear problem, at the final time t final = 27r, 
for the values e = 1 and e = 0.01. On the top line of Figure 11, we represent the plot of 
the reference solution /^^^ computed with the splitting scheme. Note that, for e = 0.01 
(top-right plot of the figure), the solution cannot be distinguished from the solution of 
the limit model, which is simply the initial data rotated of an angle 7r/2 (compare with 
Figure 1). Indeed, if oj is given by (2.35), one has ujtfinai = f + gl^ ~ f • 

On the middle line of the same figure, we represent the 2D plot of the difference 
/^p — /^g£, where /^p is the numerical solution with the scheme AP with correction 
(N = 256), for e = 1 and e = 0.01 and, on the bottom line of the figure, we represent 
the difference /second ~ /ref where /Jgcond analytic solution of the second order 

model, for the same values of e. We observe the following facts. The error for the AP 
scheme is almost the same (around 10~^) for the two values of e, whereas for the second 
order model, the results are very dependent of e: for e = 1, the error in L°° norm is 
close to 1, whereas for e = 0.01, this error is around 10~^. The second order model 
can be used only for small values of e (its incapacity to predict the solution for e = 1 
is even clearer below on the RMS test). 
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Figure 11: 2D plots of f^{t,r,v) for {r,v) £ [—2,2]^ for the linear beam model at time 
t = 27r. Top line: reference solutions computed with the splitting scheme, for e = 1 and 
e = 0.01. Middle line: difference between the reference solutions and the numerical 
solutions with the scheme AP with correction. Bottom line: difference between the 
reference solutions and the numerical solutions with the scheme second order model. 
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Accuracy of the limit model and of the second order model: Figure 12 

Let us confirm more quantitatively the above observations. In the next table, we give 
the relative L°° errors between the approximate solutions and the reference solution 
(still at time t final = 27r and, for the AP scheme, we take N = 256). This error is 
defined by 

ll/approx ~ /ref ll-f-"" 

error = . 



e 


1 


0.5 


0.25 


0.1 


0.01 


error for AP with correction 


1.8% 


1.5% 


1.5% 


1.4% 


1.3% 


error for the second order model 


18% 


4% 


1% 


0.15% 


0.001 % 


error for the limit model 


37% 


18% 


8.6% 


3.3% 


0.3% 



This table indicates that the error produced by the scheme AP with correction is 
independent of e (as we shown in the previous subsection for the nonlinear case), and 
that the limit model and the second order model seem respectively of orders 1 and 2 
in £. On Figure 12, we illustrate numerically the accuracies of these two asymptotic 
models with respect to e by plotting in logarithmic scales the L^{[0,t final]) norm of 
the difference RMSapproxit) — RM Srcfcrcncdt) , for these two models. One can check on 
this figure that the errors produced by these models are respectively 0{e) and O(e^). 
In other terms, we confirm numerically the results given by Proposition 2.4. 




epsilon 



Figure 12: Plot of the errors between the limit model and the reference solution and 
between the second order model and the reference solution, as functions of e. 

Evolution of the RMS: Figures 13, 14 and 15 

We now observe the evolution in time of the RMS quantity defined by (4.2). On Figures 



29 



13, 14 and 15, we represent respectively, for e = 1, e = 0.25 and e = 0.05, the time 
history of RMS{t) computed by the AP scheme in the linear case, with and without 
correction, and compare these numerical solutions to a reference solution and to the 
solutions of the limit model and of the second order model. In all these simulations, 
we take A'^ = 64 and At = 0.02. In all the cases, one can observe that the solution 
obtained by AP with correction fits very well with the reference solution (see the zooms 
on the right part of each figure) even when the oscillation is not well resolved. As for 
the nonlinear case, one observes that the solution obtained with AP without correction 
is less accurate when e is small and is not able to reproduce the details of order 0(e). 

One also observes that the limit model is only able to give the averaged behavior of 
the curve. The second order model is much better and follows the oscillations for small 
values of e. On Figure 15, for e = 0.05, its solution coincides with the reference solution 
and is more precise than AP with correction. Recall indeed that the error made by the 
scheme AP with correction is proportional to At^ + A^^ ~ 0.02 whereas the error made 
by the second order model is proportional to w 0.002. Indeed, in this linear context, 
the second order model is analytic and does not produce any error in time or space. 
Obviously, in a more general case, the second order model will also generate an error 
due to its space-time discretization. On Figure 14, for s = 0.25 (e^ ~ 0.06), the errors 
made by the two methods AP with correction and second order model are comparable. 
Finally, on Figure 13, for e = 1, it appears again that the error made by the second 
order model is of order 0(1): this confirms that this averaged model is useless when e 
is not small. 




time 



Figure 13: Time history of RMS(t) for e = 1, in the linear situation, computed with 
AP (with or without correction; these two curves red and blue coincide), with the limit 
model and with the second order model. On the right: zoom of the left figure. 
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Figure 14: Time history of RMS{t) for e = 0.25, in the hnear situation, computed 
with AP (with or without correction), with the limit model and with the second order 
model. On the right: zoom of the left figure. 




time time 

Figure 15: Time history of RMS{t) for e = 0.05, in the hnear situation, computed 
with AP (with or without correction), with the limit model and with the second order 
model. On the right: zoom of the left figure. 
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5 Conclusion 



In this work, we have presented a general strategy to construct the so-called Asymp- 
totic Preserving (AP) numerical schemes for a family of highly oscillatory problems. 
Although we focus on the particular case of a charged particle beam to illustrate our 
strategy, the approach may be applied to many known physical models belonging to 
this family of highly oscillatory problems. Averaged models are usually used to ap- 
proximate this type of problems but these models are not relevant in the intermediate 
regime since they miss important informations from the original problem. 

The starting idea in this construction is to write the oscillatory problem into a 
"double-scale" formulation where the rapid and slow time scales are separated, making 
the new distribution function more regular in some sense. The new structure then 
suggests to follow a similar strategy as in the collisional case to develop AP schemes on 
this formulation. However the completely different nature of highly oscillatory problems 
(compared to collisional kinetic equations) induces new important difficulties. First, the 
double-scale formulation is overdetermined in the sense that a large family of initial data 
for this formulation is allowed. We show in this paper that there is a suitable choice to 
make on this initial data in order to maintain the regularity of the distribution function 
at different orders of the oscillation parameter. More precisely, the initial data is chosen 
to fit with a Chapman-Enskog like expansion which ensures a separation of the rapid 
and slow time scales at different orders of the expansion. Based on this formulation, 
we then derive an Asymptotic Preserving scheme for the original problem and show 
that time-space discretizations of order 2 are necessary to numerically observe the fine 
structures and filamentations that are generated by the coupling of Vlasov and Poisson 
equations. Several numerical tests are performed to show the efficiency of our strategy: 
uniform accuracy and ability to capture the oscillations of different magnitudes and the 
long time behavior. 

We emphasize that the AP property of our scheme is shown by making links with the 
so-called micro-macro decomposition, which is known to be a flexible tool to develop AP 
schemes in the context of collisional kinetic equations. In particular, this decomposition 
may be used to extend the present approach to other highly oscillatory problems such 
as the charged particle beam with a diffusion scaling, the guiding-center asymptotics 
and the finite Larmor radius approximation. This will be the subject of future works. 

References 

[1] G. Allaire, Homogenization and two-scale convergence, SIAM J. Math. Anal. 23, 
1482-1518, 1992. 

[2] M. Bennoune, M. Lemou, L. Mieussens, Uniformly stable numerical schemes 
for the Boltzmann equation preserving compressible Navier-Stokes asymptotics, J. 
Comput. Phys. 227, 3781-3803, 2008. 

[3] N. Besse, E. Sonnendrucker, Semi-Lagrangian schemes for the two-dimensional 
Vlasov-Poisson system on unstructured meshes, J. Comput. Phys. 191, 341-376, 
(2003). 



32 



[4] M. BOSTAN, The Vlasov-Poisson system with strong external magnetic field. Finite 
Larmor radius regime, Asymptot. Anal. 61, 91-123, 2009. 

[5] M. BoSTAN, The Vlasov-Maxwell system with strong initial magnetic field. Guiding- 
center approximation, Multiscale Model. Simul. 6, 1026-1058, 2007. 

[6] F. Castella, p. Chartier, A. Murua, F. Mehats, Stroboscopic averaging for 
the nonlinear Schrodinger equation, preprint HAL-00732850 (http://hal.archives- 
ouvertes.fr). 

[7] N. Crouseilles, M. Lemou, An asymptotic preserving scheme based on a micro- 
macro decomposition for collisional Vlasov equations: diffusion and high-field scal- 
ing limits. Kinetic Related Models 4, 441-477, 2011. 

[8] N. Crouseilles, E. Frenod, S. Hirstoaga, A. Mouton, Two-Scale Macro- 
Micro decomposition of the Vlasov equation with a strong magnetic field, to appear 
in Math. Models Methods Appl. Sci. 

[9] P. Degond, p. a. Raviart, The paraxial approximation of the Vlasov-Maxwell 
equations, Math. Models Methods Appl. Sci. 3, 513-562, 1993. 

[10] F. Filbet, E. Sonnendrucker, Modeling and numerical simulation of space 
charge dominated beams in the paraxial approximation. Math. Models Methods 
Appl. Sci. 16, 763-791, 2005. 

[11] E. Frenod, P. -A. Raviart, E. Sonnendrugker, Two scale expansion of a 
singularly perturbed convection equation, J. Maths. Pures Appl. 80, 815-843, 2001. 

[12] E. Frenod, F. Salvarini, E. Sonnendrucker, Long time simulation of a beam 
in a periodic focusing channel via a two-scale PIC-method, Math. Models Methods 
Appl. Sci. 19, 175-197, 2009. 

[13] E. Frenod, E. Sonnendrugker, Long time behavior of the Vlasov equation with 
strong external magnetic field. Math. Models Methods Appl. Sci. 10, 539-553, 2000. 

[14] E. Frenod, E. Sonnendrugker, The finite Larmor radius approximation, SIAM 
J. Math. Anal. 32, 1227-1247, 2001. 

[15] F. GOLSE, L. Saint-Raymond, The Vlasov-Poisson system with strong magnetic 
field, J. Math. Pures Appl. 78, no. 8, 791-817, 1999. 

[16] A. R. GOURLAY, J. L. Morris, Multistep formulation of the optimized Lax- 
Wendroff method for nonlinear hyperbolic systems in two space variables. Math. 
Comp. 22, 715-719, 1968. 

[17] S. Jin, Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic 
equations, SIAM J. Sci. Comput. 21, 441-454, 1999. 

[18] M. Lemou, Relaxed micro-macro schemes for kinetic equations, C. R. Math. Acad. 
Sci. Paris 348, 455-460, 2010. 



33 



[19] M. Lemou, F. Mehats, Micro-macro schemes for kinetic equations including 
boundary layers, to appear in SIAM J. Sci. Comput. 

[20] M. Lemou, L. Mieussens, A new asymptotic preserving scheme based on micro- 
macro formulation for linear kinetic equations in the diffusion limit, SIAM J. Sci. 
Comp. 31, 334-368, 2008. 

[21] J.-G. Liu, L. Mieussens, Analysis of an asymptotic preserving scheme for linear 
kinetic equations in the diffusion limit, SIAM J. Numer. Anal. 48, 1474-1491, 2010. 

[22] A. MOUTON, Two-scale semi-Lagrangian simulation of a charged particles beam 
in a periodic focusing channel, Kinetic Related Models 2, 251-274, 2008. 

[23] G. Nguetseng, a general convergence result for a functional related to the theory 
of homogenization, SIAM J. Math. Anal. 20, 608-623, 1989. 

[24] L. M. Perko, Higher order averaging and related methods for perturbed periodic 
and quasi-periodic systems, SIAM J. Applied. Math. 17, 698-724, 1969. 

[25] R. D. Richtmyer, A survey of difference methods for non-steady fluid dynamics, 
NCAR technical notes, 1962. 

[26] J. A. Sanders, F. Verhulst, Averaging methods in nonlinear dynamical sys- 
tems. Applied Mathematical Sciences, Vol. 59. Springer- Verlag, 1985. 

[27] X. Yang, F. Golse, Z. Y. Huang, S. Jin, Numerical study of a domain decom- 
position method for a two-scale linear transport equation. Networks and Heteroge- 
neous Media 1, no. 1, 143-166, 2006. 



34 



